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Abstract. A class of parametric distribution functions has been proposed in 
[C. DiTroia, Plasma Physics and Controlled Fusion, 54, (2012)] as equilibrium 
distribution functions (EDFs) for charged particles in fusion plasmas, representing 
supra-thermal particles in anisotropic equilibria for Neutral Beam Injection, Ion 
Cyclotron Heating scenarios. Moreover, the EDFs can also represent nearly isotropic 
equilibria for Slowing-Down alpha particles and core thermal plasma populations. 
These EDFs depend on constants of motion (COMs). Assuming an axisymmetric 
system with no equilibrium electric field, the EDF depends on the toroidal canonical 
momentum the kinetic energy w and the magnetic moment p. 

In the present work, the EDFs are obtained from first principles and general hypothesis. 
The derivation is probabilistic and makes use of the Bayes’ Theorem. The bayesian 
argument allows us to describe how far from the prior probability distribution function 
(pdf), e.g. Maxwellian, the plasma is, based on the information obtained from magnetic 
moment and GC velocity pdf. 

Once the general functional form of the EDF has been settled, it is shown how to 
associate a Landau collision operator and a Fokker-Planck equation that ensures the 
system relaxation towards the proposed EDF. 


PACS numbers: 52.55.Fa, 52.65.Ff, 02.50.Cw 
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1. Introduction 

The following parametric distribution function has been proposed in Ref.[T] as 
equilibrium distribution function (EDF) for charged particles in fusion plasmas, 
representing, e.g., supra-thermal particle distribution produced by additional external 
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heating sources in tokamak experiments: 
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being w the kinetic energy per unit mass, /i the magnetic moment per unit mass, 
A = /i/tc the pitch angle and = {es/ms)P(f>, being e* the charge and nis the mass of the 
considered species, and p^ the canonical toroidal momentum, assuming an axisymmetric 
system. Moreover, AA, T^„, Ap_^, Aq and are control parameters. In pQ, the 

orbit theory has been described through the constant of motions (COMs) w, A, where 
the canonical momentum is treated as a spatial coordinate; the same choice is taken 
also her^. 

Together with ([1]), the regularized EDF, 

feq,R feqheqi (2) 

being 

11(1^^ Wj)( 5 confined / q \ 

has been proposed in [1] as general plasma EDF for describing populations of particles 
encountered in many tokamak scenarios. In ([3]), H(wfe —w) is the Heaviside step function 
which takes into account the presence of a mono-chromatic source of birth energy Wb. 
The factor 1 -f {wc/wY^"^ mimics the Slowing-Down behavior in energy, being Wc the 
critical energy 0 [in], resulting from the relaxation of the considered species with bulk 
ions and electrons. The symbol (^confined is the analytical condition for a particle whose 
orbit is mostly determined by and A, to be conhned in the plasma volum^. 

This EDF has already been implemented in the hybrid code XHMGC |2] and in the 
gyrokinetic code NEMORB [3]. It has been shown that, by varying the EDF control 
parameters AA, au,,T^, Ap^, Aq and A;i,, (|2]) can represent anisotropic equilibria as for 
the case of Neutral Beam Injection and Ion Cyclotron (or Electron Cyclotron) Resonance 
Heating. Moreover, it can also represent nearly isotropic equilibria as for the case of 
Slowing-Down alpha particles and core thermal plasma populations. In jl] it has been 
proposed a heuristic derivation of feq whilst, in the present work, a rigorous one is 
shown based on probabilistic principles and general hypothesis for deriving a class of 
EDFs which includes also the distribution function ([T]) and ([ 2 ]). 

The distribution function, /eg = /eg(P 0 , w, A), in ([1]) is an EDF because it depends solely 
on COMs. In this way the total time derivative is 

feq = V<t>9p^feq + wd^feq + XOxfeq = 0, (4) 

being P^ = re = A = 0. In an EDF, the dependency on COMs is commonly obtained 
by the transport Boltzmann equation, and precisely by the kernel of the Boltzmann 
collision operator, Cp: 

CBUeq) = 0. (5) 


f At equilibrium the motion is unpertubed and fields are stationary, so that it will be considered only 
the guiding center transformation. 

§ The explicit analytical expression of (^confined can be found in [1] and will not be reported here. 
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Often the problem for solving the equilibrium of the Boltzmann equation for plasmas is 
attacked in the following manner: 

1) give an analytical expression to the collision operator, 

2) hnd the distribution function that belongs to the kernel of the collision operator, 

3) check if a combination of such solutions is constant in time and 

4) try to express such combination of solutions as the function of a particular set of 
COMs. 

A simplihcation of the problem is realized just starting from the maxwellian distribu¬ 
tion function which is the known solution of the Boltzmann equation, expressed with 
the Landau collision operator, for a gas of charged particles interacting via Coulomb 
collisions. The maxwellian EDF is further transformed into a local maxwellian and, 
further, into a canonical maxwellian to end on writing the EDF as a unbiased canonical 
maxwellian |1]. 

The problem is that the maxwellian satishes the Boltzmann equilibrium for Coulomb 
interactions but it is not function of COMs, energy apart. Moreover, although the unbi¬ 
ased canonical maxwellian is constant in time, it doesn’t really belong to the standard 
Landau collision operator kernel. 

Here it is proposed an alternative way to reach an EDF, using the Bayes' theorem; the 
maxwellian distribution function is considered only as the prior probability distribution 
function (pdf), while the joint distribution function to have a certain probability to hnd 
a particle with given w and A are considered the hnal EDF expressed as the product 
of the conditional probability, to have such Vij, and A once the energy of the particle 
is known, multiplied for the prior pdf. The obtained EDF, being different from the 
maxwellian EDF, cannot be anymore considered as the solution of the standard Landau 
collision operator but it will be shown that it belongs to the kernel of a Boltzmann 
collision operator which is a modihed Landau collision operator. 

The problem of Ending an EDF depending entirely and solely on COMs is addressed 
in section 2, while the form of the associated Landau collision operator is studied in 
section 3. 


2. Probabilistic derivation of the equilibrium distribution function 

The position a: of a charged particle in a magnetic held is expressed by 

X = X p, (6) 

where X is the Guiding Center (GC) and p is the vector Larmor radius. Similarly the 
particle velocity v is 

V = X = V + a, (7) 

where D = X is the GC velocity and a = p is the diherence between the particle and 
the GC velocities. In the GC transformation, the velocity is expressed hy v = v\\b + v 
separating a component parallel to the magnetic held B [b is the unit vector of B) to 
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the component perpendicular to it. At first it is possible to associate V with v\\b and a 
with ux, but this correspondence can be done only approximatively, being exactly only 
for constant and uniform magnetic held when the drift velocity is zero. In general V 
and a are not orthogonal vectors and the angle a between them is important for the 
present derivation: 


cos a = 


a-V 

din’ 


( 8 ) 


being |l/|s = sgn(l/ ■ B)\V\, where the sign of V depends on its orientation towards B. 
The angle a belongs to [0, tt] if V y' B, and a G [tt, 27t] if V \ B. |1/|s generalizes the 
parallel velocity uy, and it is particularly useful when the GC transformation is taken 
at all perturbative orders as done in [5]. 

The value cos a is considered as the realization of the random variable Cosa. Such 
assumption can be also considered a consequence of the GC transformation which treats 
the gyro-angle 7 as an ignorable coordinate as if it is possible to randomly chose a value 
of 7 without affecting the equations of motion in GC coordinates. This can be done 
also for cosa because it depends strictly on 7 . The gyro-angle 7 can be assumed to 
be uniformly distributed between [0, 27r), whilst the distribution of cosa will depend on 
the considered physical scenario. Regardless of the scenario, it is assumed that such 
distribution is centered at zero. A simple situation where the value of cosa is always 
zero, is for constant and uniform B, when V is orthogonal to a. 

It is here assumed that Cosa is distributed as a Normal distribution function with a 
deviation depending on the parameter k. In this way, given all particles with same GC 
position X, same value of the GC velocity |R|s and same kinetic energy w, then 


Cosa ~ N(0,k/2|1/|), 


( 9 ) 


where k will be determined later, at the end of the present section, to discriminate 
various tokamak scenarios. 

From the knowledge of w, |R|s and cosa, it is possible to compute |(t| as below 
shown starting from the equivalence 2w = = (V + a)^. IndeeclJJ], 

-|-2|ct| |R|s cos a — (2t(; — = 0, (10) 


with the solution 


a\ = —\V\s cosa -|- cos^ a + 2w — V^, 


( 11 ) 


rewritten for convenience as 
|(t| = \/2w — V^ 


|R|sCosa 
^/2w - R2 


cos^ a 
2w-V^ 


+ 1 


( 12 ) 


As before, |(t| can be considered the realization of a random variable because it is 
function of cosa. The same can also be said for the following variable 

m = -, 

2msUc 



II The present derivation can also be used on describing elementary particle decaying processes if the 
products are a pair of particles departing at right angles. 
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where cOc is the cyclotron frequency, dehned as Uc = 7 . The variable m is an estimate 
of the magnetic moment /r, here dehned as: 
w-V^l2 

/i = - . (14) 

msOJc/(^s 

The choice on the above, unusual dehnition of the magnetic moment, comes from having 
considered ignorable the gyro-phase 7 and taken a null electric potential. In such case 
the single particle hamiltonian is the kinetic energy expressed in the canonical variables: 
(X, P = V + CsA/ma) and {cs'j/ms, A being the vector potential, so that 

(P - CsAImff 


w = 


+ p{msOJc/e.s), 


(15) 


from where the dehnition (ITTl) is taken. It is worth noticing that, from the above 
dehnition of the magnetic moment, p is an exact COM. Indeed, the Hamilton's equations 
are: {es/ms)^ = d^w = Uc and p = {es/ms)d^w = 00 

The variable m is a good estimate of p when cos a ~ 0. From equation flT^ . the random 
variable M, of which the value m is a realization, is explicitly written as 

2 


M = p 


kY 

-L . 

( )''' 

2V2w -V^ \ 

\2V2w-Vy 


being 


Y = 


2\V\sCosa 
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(16) 


(17) 


From ([9]), the random variable Y is distributed as Y N(0,1). 
Introducing a and (3\ 

..2 
, 2 


l3 = p and a = 


K 


2p{msUJcl Cs) 


then M is rewritten as 


M = 13 



(18) 


(19) 


which is recognized as the Birnbaum-Saunders (BS) [ 6 ] random variable: M ~ BS(a;, /3). 
The corresponding pdf, also known as fatigue life pdf, is 


fBs{m;a,P) = 


(3 + m 


exp 



( 20 ) 


(27r/3)^/2 2aim?P 

The BS pdf has been developed to model breakage due to cracks, when a material is 
subjected to cyclic stress and every Tth cycle leads to an increase in crack extension of 
Y = Yi. This crack grows under the repeated applications of a common cyclic stress 
pattern, until it reaches a critical size, when fatigue failure occurs. If the total extension 
of the crack is normally distributed, then T, the time until failure, is distributed as a 

^ When the electric potential, $(t,a;), must be considered, it occurs to substitute the kinetic energy 
with the total particle energy: w ^ e = w + {es/ms)^{t,x) = F^/2 + {es/ms)^{t, X) + y{msUJc/es), 
as described in [5]. 
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The pdf of m is a conditional one, because a and fd are supposed to be already 
known, depending on w, k and X. A similar argument can be applied for the 

estimate V of |l/|s. Suppose to be able to make a measurement, affected by an error, of 
the GC velocity of a particle and repeat such measurement for all particles with same 
iGls and same GC position X. The result V is assumed to follow a Gaussian behavior 
with a standard deviation Ay/\/2: 


fG{V-\Vs\,Av/V2) 


1 



(^) 


( 21 ) 


The system is considered in the Maxwellian state as prior pdf, expressed by the following 
exponential law: 


gM{w]Tyf) = — exp 

J- 'in 


w 

T~ 


( 22 ) 


This assumption is standard, but can also be generalized considering the Gamma 
distribution: 


^p('U), CT-U,, Tyfj 


{w/Tu 


\ 1 


exp 


w 


(23) 


nr{a^ 

Some inferences are required to know how far the system is from being described by the 
prior distribution. The desired pdf is obtained from the conditional probability 7r(6* | x) 
to have 9 once x is given, where 9 is the array of true values 9 = {\V\s,w,fi) and x the 
array of the “estimated” quantities x = {V,m). From Bayes'' theorem: 

TT{9)f{x I 9) 


n{9 I x) = 


(24) 


I 6')dv(e')’ 

which means that the posterior probability, 7t{ 9 \ x), is proportional to the product of 
the prior n{9), fl2^ or fl23l) . multiplied for the conditional probability fl20|) and fl2l]) : 

A(/i + m) 


■k{9 I m, G) = 


X 


Tiy/^ m^/^aTy^Av 

( 1 ^ 1 . 


The above distribution function is not a function of COMs, but it is possible to properly 
re-cast it into an explicit EDF form, hrstly, by substituting the space parameters 
0 —)■ ©COM = (7^4), w,p,). This transformation can be realized from the linear 


(/i — mf 


+ For plasmas, the analogy is interesting: the width of the crack is associated to the weighted projection, 
depending on the scenario, of V on tr. Such extension is zero for constant and uniform magnetic fields. 
The random variable, T, is assumed distributed as a Normal pdf, and one particle corresponds to one 
cycle as well as one crack extension, while the total extension corresponds to the sum of Y over the 
particles with assigned |F|s and w. The probability density function of m, whose mean value is the 
magnetic moment, follows the same behavior of the time until failure of the stressed material. 
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dependency of the difference and the magnitude of GC velocity divided for 

up. 

.\V\s 


= 7F- 


Un 


(26) 


The function of proportionality, can be computed explicitly from the GC 
transformation [5]. At lowest order, 7F is F, dehned in i? = x V0 + FV0. Indeed, at 
this order, |G|s —)■ ny and — = Fv\\/uc. The constant parameter V^po is conveniently 

dehned to be 


1/ 

'PfpO = V' + ^—• 

Uc 

The pdf becomes 

7r(P0,tc,/i I m, G) = 


A^(p + m) 


7ri/2p 
X exp 
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The above distribution function represent an equilibrium only if also a and m are COMs. 
In terms of v? and m the following cases are considered. 

As hrst choice. 


2 A2(m,/e,)a;c , 

K = — - - and m = So 


m 


so that 


becomes 


feql{V^,W,n) = N X 


(29) 
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This case can be useful when the system is artihcially prepared to select only particles 
with a magnetic moment sq or, better, for a classic gas of quantum charged particles 
with spin proportional to sq. If the spread on the magnetic moment is sufficiently 
large, it can reasonably represent the bulk populations of plasma at thermal equilibrium. 
This simple distribution function seems to be original as not yet proposed in literature. 
As second choice. 


2 T^^\{ms/es)upw 
K = - and m = Aotc 


m 


so that 
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^ 7^ ... ..N N{\o+fv/w){w/T^) 
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(31) 
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Substituting X = p./w, the distribution function is almost identical with ([T]) apart from 
the multiplicative factor (1 + A/Aq), which is of minor importance in comparison to 
the exponential behavior. It is worth noticing respect to [1] that here the derivation is 
probabilistic. Moreover, the constants of motion are exacts (in ([1]) they were computed 


0^7, 


= 1 . 


only at leading order) and the power factor, in ([T]) is a direct result; 

Obviously, if another prior pdf is used, the hnal EDF will change; e.g. if (12311 is used 
instead of fl2^ . then the EDF becomes general as CO. being a^, no longer hxed. It is 
worth noticing that if a Slowing Down distribution function is used as the prior, being 

TsSa H(m;6 


fsD{w) = 


W) 


(33) 


8y/27r Wc^'^ 

where ts is the Spitzer SD time [S], is the birth energy, Wc = vl/2 is the critical 


energy, and Sa is a normalization constant, then, from 

^TsSa{l + fi/so) II(wb - w) 


feq^{V^,W,ti) = N- 


8'/27i‘^AnAp + Wc^‘^ 


the EDF becomes 


X 


(34) 


X exp 
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(h - so)^ 
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Otherwise, from (|3Tll . it is obtained: 

.tsS^{1 + X/\o){w/T^)-^/^ H(n;fe 
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(35) 


already considered in [T] for representing fast particles heated by Neutral Beam Injec¬ 
tion. 

Other EDFs can be obtained by varying some assumptions made here concerning the 
random variable Y. Indeed, the BS is a particular pdf belonging to the family of Inverse 
Gaussian pdfs [HI HSl [13] and it can be generalized as already studied in [HI [15]. 

In the present derivation it is also furnished the method for consistently addressing the 
values of the control parameters. The parameters Aq and V,po are obtained from the 
estimate of the magnetic moment and the estimate of the magnitude of GO velocity, m 
and V, respectively, which characterized some particular orbits thus becoming represen¬ 
tative of the considered population. If such quantities cannot be measured or estimated, 
it is not a problem because the functional form of the EDF is maintained. In this case, 
the EDF will be mainly used as a htting model function. This doesn’t mean that the 
parameters are devoid of any physical interpretation. It simply means that such param¬ 
eters should be inferred from what is measurable as, for example, density, temperature 
or pressure of the considered plasma species might be. 
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The Boltzmann equation for a distribution function fs = /^(t, x, n), in which s indicates 
the species of particles with mass and charge e^, is fg = CB{fs), where Cb is the 
Boltzmann collision operator. An important class of these operators consists of the 
following Fokker Planck operators: 


Cppifs) — • {Dg ■ Vv + dg)fg, 


(36) 


being Dg the diffusion matrix and dg the collisional drag. In such case, the Boltzmann 
equation is said Fokker-Planck equation and applies for describing soft collisions, i.e. 
binary collisions with only little transfers of velocity changes to scattered particles. If a 
background species of mass m^/ and charge e^', is recognized to work as scatterers with a 
distribution function /^, = /(t, x', n'), then the scattering Fokker-Planck operator would 
be Cpp^fg) = J2s' Cppifgi, fg) where the diffusion matrix and the collisional drag are 
functional of f^r. Dg = Dg{fg,) and dg = ds{f'^,). A smart representation of the above 
operator is given by the Landau collision operator: 

CMpfs) = • / d^v'U{u) • - fgVpfg,) . (37) 

being '^g/g a constant which is known for the Coulomb collision case, and where 

Jd\'U(u)f:„ (38) 

and 


4..(/;) = ^/rfVr/(«)-v„./),. 

The scattering matrix U is 



(39) 

(40) 


for u = u{v, v') that will be specihed later. A useful features of Landau collision operator 
is that it can be described by the RMJ potentials [16]. Moreover, the main characteristic 
of Cl is the assurance of the relaxation system to equilibrium thanks to the entropy 
production, generally dehned by 


Q = - J d^vCp log fg ; 


(41) 


it can be shown that 0 > 0 for any fg which satishes the Boltzmann equation written 
with the Landau collision operator. Moreover, if e.g. Cp is adopted, the equilibrium, 
/eg, is reached only if Cpifeq) = 0 when 0 = 0. 

The standard application of Cl for deriving the EDF for a dilute plasma via 
Coulomb interactions leads to the maxwellian distribution function, /m: 


/m — 


n 


exp I —— 


w 

T 

-L 11 


(27rT^)3/2 

Indeed, if fM and are substituted in fl37ll for same species and same T^, then 


(42) 


CL{f'MJM) = 


7c 

2T 

in 


d^v' U(v, v') ■ (v - ii')/m/m, 


(43) 
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being jc constant for Coulomb collisions. Now, if m = n — n' is the relative velocity 
between two colliding particles then U ■ {v — v') =0 and also Cl = 0. Such result 
led many people to the belief that the maxwellian distribution function is the only 
equilibrium distribution function allowed for a plasma. 

Hence, a simple method to overcome this common belief and allows the EDF fl30l) 
or ([1]) to represent a plasma in a magnetic held at equilibrium is now considered. The 
problem is resolved appropriately changing the Landau collision operator. The velocity 
vector M, taken as the relative velocity for the maxwellian equilibrium, changes to Upi, 
dehned as 



(44) 


and leaving the same U matrix in the same Cl operator. From (l44ll . 



\_Uc /i + So w' /i' + So 


at same temperature, = T^, same magnetic moment spread, A^, and same mean 
magnetic moment, So. The same Landau collision operator and the corresponding 
entropy production, bring the system to the EDF fl3UI) . which is maxwellian-like when 
A^ and so are sufficiently large. The maxwellian limiting behavior is recognized also 
in (05]). Indeed, if A^ and s are large, then Upi comes to be the same relative velocity 
which brings the system to the maxwellian EDF. It is also worth noticing that the 
high Uc frequency at the denominator in fj45ll implies that the equilibrium is close to 
a maxwellian also when the magnetic held is sufficiently strong or when the Larmor 
radius, pl = |a|/|a;c|, is sufficiently small. Although such model doesn’t accurately 
to describe the true interactions between charged particles, it can be suggested that 
the velocity vector Upi in (I44|) is the sum of the relative velocity, v — v', that takes into 
account Coulomb collisions, and a new part representative of other interactions, as it can 
be the Ampere''s interactions between currents (indeed, both a as a' are the velocities 
of charged particles moving on closed loops). Moreover, if a —?■ 0 then Up = V — V\ 
becomes the relative velocity of colliding guiding particles, that are particles with the 
same velocity and position of GC, and null Larmor radius [5]. This means that the 
equilibrium distribution of guiding particles is, again, maxwellian with energy V‘^/2. 
Same considerations are true if the EDF in ([1]) is chosen to represent the thermal 
population equilibrium with the only difference to compute U ioi u = Up\ 



(46) 



2 


2 (A^ - Ag) - AjT^/n; ] 

A + Ao J 

' 2 (A-^ - Ag) - A.lT^/w' - 
A' +Ao 
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with same EDF’s parameters, T^, Aq, and (the latter parameter derives from the 
prior distribution in fl23|l instead of the fl2^ h The above expression is more complicated 
than (05]), which is the preferred choice, although in this context, it is not possible to 
exclude (06|) . It is worth noticing that if —)■ oo and = 3/2, then 


Up ^ V — v' — (a — a’'' 


e.,r„ 


msWUJc{A + Aoj 

and it becomes, again, easy to represent a maxwellian-like EDF, for {mg/es)u;c{^ + Xq) 3> 
(T^/w). 

Unfortunately, a choice is necessary because the same stratagem of modifying the 
Landau collision operator cannot be applied twice for deriving two or more EDFs. 
However, there is a difference between the EDF in fl30D (or in ([T|)) and the others, 
dlj), fIMD and (I5BD . Indeed, the lasts are describing plasma populations in the presence 
of sources and losses (of particles or energies). The equilibrium distribution function 
feq,s G {feq,R, feq3, feq4:}, IS obtained from the balance between the collisions with the 
background populations, i. e. bulk ions and electrons represented by ft, G {feqi (or 
feq}), and the source and the losses of particles and energies: 


0 = feq,s = CB{feq,s) = E ^dfi feq,s) + S, + L,. (48) 

b 

The species labelled with s can be alpha products. He minority, energetic Deuterium, 
etc. In this way, the losses and the sources are modeled with the same parameters of 
the EDFs; indeed, 

S, + L, = - E ^ V, • J d\' U{u) ■ {nvjeq,s - feq,sV,,Q , (49) 

b 

where the EDF of the bulk is = /4i(or fl, = f^q), and u = Mpi(or u = Up). This 
remark is quite important for tokamak plasmas because some of the sources, e.g. radio¬ 
frequency antennas or neutral beam injectors, are controlled and can be appropriately 
modeled to be consistent with the desired bulk equilibrium. 


4. Conclusions 

In the present work, the probabilistic derivation of axisymmetric plasma EDFs, making 
use of the Bayes^ theorem, has been explained. Four EDFs, fl30l) . fl32ll . fl34l) and fl35ll . are 
obtained explicitly depending on if the mean magnetic moment is almost constant for 
particles with different energies, or if it scales in magnitude with the particle energy, and 
if a Maxwellian or a Slowing Down is considered as prior pdf. The pdf fl5UD can represent 
a plasma in thermal equilibrium in a particular limit behavior, when the parameters A^ 
and So are sufficiently large. The pdfs fl32l) and fl35l) are recognized to be the same EDF 
already proposed in Ref. [1] which is known to be useful when external heating sources, 
i. e. Ion Cyclotron Resonance Heating or Neutral Beam Injection, are employed ini 
in experimental tokamak campaigns. The pdf (j3T|l represents a more realistic EDF for 
fusion products in axisymmetric tokamak. 
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The relevance of the present derivation resides in its generality since it only reqnires the 
GC transformation of phase space coordinates and an axisymmetric ambient magnetic 
held. Moreover, snch derivation doesn’t depend on the detailed form of the axisymmetric 
magnetic held. This means that the derived fnnctional form of the EDF is machine 
independent. Only the COMs and /i are fnnctions of the eqnilibrinm B so that 
diherent valnes of the EDF parameters correspond to diherent scenarios. Respect to 
|1], where only the leading order approximation of COMs in GC coordinates are nsed, 
in this case instead, the considered COMs, V^, w and ja are all exact invariants. In 
conclnsion, the proposed EDFs is placed in a more ’’thermodynamic” framework with 
a little change of the Landau collision operator, which is still preserved in the form. 
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